function eigenvalue_basics()
% EIGENVALUE_BASICS 特征值和特征向量基础
% 
% 介绍特征值和特征向量的基本概念、性质和计算

fprintf('=== 特征值和特征向量基础 ===\n\n');

%% 1. 基本定义
fprintf('1. 基本定义\n');
fprintf('对于n×n矩阵A，如果存在标量λ和非零向量v使得 Av = λv\n');
fprintf('则λ称为特征值，v称为对应的特征向量\n\n');

% 简单示例
A = [3, 1; 0, 2];
fprintf('矩阵 A =\n');
disp(A);

% 计算特征值和特征向量
[V, D] = eig(A);
eigenvals = diag(D);

fprintf('特征值: λ₁ = %.6f, λ₂ = %.6f\n', eigenvals);
fprintf('特征向量矩阵 V =\n');
disp(V);

% 验证特征值方程
for i = 1:length(eigenvals)
    lambda = eigenvals(i);
    v = V(:, i);
    Av = A * v;
    lambda_v = lambda * v;
    
    fprintf('验证 λ_%d: ||Av - λv|| = %e\n', i, norm(Av - lambda_v));
end

%% 2. 特征多项式
fprintf('\n2. 特征多项式\n');
fprintf('特征多项式: det(A - λI) = 0\n\n');

% 符号计算特征多项式（简化版本）
fprintf('对于2×2矩阵 A = [a, b; c, d]\n');
fprintf('特征多项式: λ² - (a+d)λ + (ad-bc) = 0\n');
fprintf('即: λ² - trace(A)λ + det(A) = 0\n\n');

trace_A = trace(A);
det_A = det(A);
fprintf('trace(A) = %.6f\n', trace_A);
fprintf('det(A) = %.6f\n', det_A);

% 验证特征多项式的根
fprintf('特征多项式的根:\n');
coeffs = [1, -trace_A, det_A];
roots_poly = roots(coeffs);
fprintf('根1: %.6f, 根2: %.6f\n', roots_poly);
fprintf('与特征值的差异: %e\n', norm(sort(roots_poly) - sort(eigenvals)));

%% 3. 特征值的基本性质
fprintf('\n3. 特征值的基本性质\n');

% 性质1: 特征值的和等于迹
sum_eigenvals = sum(eigenvals);
fprintf('特征值之和: %.6f\n', sum_eigenvals);
fprintf('矩阵的迹: %.6f\n', trace_A);
fprintf('差异: %e\n', abs(sum_eigenvals - trace_A));

% 性质2: 特征值的积等于行列式
prod_eigenvals = prod(eigenvals);
fprintf('特征值之积: %.6f\n', prod_eigenvals);
fprintf('矩阵行列式: %.6f\n', det_A);
fprintf('差异: %e\n', abs(prod_eigenvals - det_A));

% 性质3: 相似矩阵有相同的特征值
P = [1, 2; 1, 3];
A_similar = P \ A * P;
eigenvals_similar = eig(A_similar);

fprintf('相似矩阵的特征值: [%.6f, %.6f]\n', eigenvals_similar);
fprintf('特征值差异: %e\n', norm(sort(eigenvals) - sort(eigenvals_similar)));

%% 4. 复特征值
fprintf('\n4. 复特征值\n');
fprintf('实矩阵可能有复特征值\n\n');

% 旋转矩阵示例
theta = pi/4;
R = [cos(theta), -sin(theta); sin(theta), cos(theta)];
fprintf('旋转矩阵 R (θ = π/4) =\n');
disp(R);

[V_R, D_R] = eig(R);
eigenvals_R = diag(D_R);
fprintf('特征值: λ₁ = %.4f + %.4fi, λ₂ = %.4f + %.4fi\n', ...
        real(eigenvals_R(1)), imag(eigenvals_R(1)), ...
        real(eigenvals_R(2)), imag(eigenvals_R(2)));

% 验证 |λ| = 1 (旋转矩阵的特征值在单位圆上)
eigenval_magnitudes = abs(eigenvals_R);
fprintf('特征值的模: |λ₁| = %.6f, |λ₂| = %.6f\n', eigenval_magnitudes);

%% 5. 对称矩阵的特征值
fprintf('\n5. 对称矩阵的特征值\n');
fprintf('对称矩阵的特征值都是实数，特征向量相互正交\n\n');

% 对称矩阵
S = [2, 1; 1, 3];
fprintf('对称矩阵 S =\n');
disp(S);

[V_S, D_S] = eig(S);
eigenvals_S = diag(D_S);

fprintf('特征值: λ₁ = %.6f, λ₂ = %.6f\n', eigenvals_S);
fprintf('特征值都是实数: %s\n', mat2str(all(imag(eigenvals_S) == 0)));

% 验证特征向量的正交性
v1 = V_S(:, 1);
v2 = V_S(:, 2);
dot_product = v1' * v2;
fprintf('特征向量内积: %.6e (应该为0)\n', dot_product);

% 验证特征向量矩阵是正交的
orthogonality_error = norm(V_S' * V_S - eye(2), 'fro');
fprintf('V^T*V - I 的误差: %e\n', orthogonality_error);

%% 6. 特征值的几何重数和代数重数
fprintf('\n6. 特征值的几何重数和代数重数\n');

% 重特征值矩阵
A_mult = [2, 1, 0; 0, 2, 1; 0, 0, 2];
fprintf('矩阵 A =\n');
disp(A_mult);

[V_mult, D_mult] = eig(A_mult);
eigenvals_mult = diag(D_mult);
fprintf('特征值: [%.6f, %.6f, %.6f]\n', eigenvals_mult);

% 分析重数
unique_eigenvals = unique(round(eigenvals_mult, 10));
fprintf('不同的特征值: %.6f\n', unique_eigenvals);

for lambda = unique_eigenvals'
    % 代数重数
    algebraic_mult = sum(abs(eigenvals_mult - lambda) < 1e-10);
    
    % 几何重数 = dim(null(A - λI))
    geometric_mult = size(A_mult, 1) - rank(A_mult - lambda * eye(size(A_mult)));
    
    fprintf('λ = %.6f: 代数重数 = %d, 几何重数 = %d\n', ...
            lambda, algebraic_mult, geometric_mult);
end

%% 7. 幂法预览
fprintf('\n7. 幂法预览\n');
fprintf('迭代计算主特征值和特征向量\n\n');

% 使用幂法计算最大特征值
A_power = [4, 1; 1, 3];
x0 = [1; 1];  % 初始向量
max_iter = 10;

fprintf('幂法迭代:\n');
x = x0;
for k = 1:max_iter
    x_new = A_power * x;
    lambda_approx = x' * x_new / (x' * x);  % Rayleigh商
    x = x_new / norm(x_new);  % 归一化
    
    if k <= 5 || k == max_iter
        fprintf('迭代 %2d: λ ≈ %.6f\n', k, lambda_approx);
    end
end

% 与精确值比较
eigenvals_exact = eig(A_power);
max_eigenval = max(eigenvals_exact);
fprintf('精确最大特征值: %.6f\n', max_eigenval);
fprintf('幂法误差: %e\n', abs(lambda_approx - max_eigenval));

%% 8. 条件数和特征值敏感性
fprintf('\n8. 条件数和特征值敏感性\n');
fprintf('特征值对矩阵扰动的敏感性\n\n');

% 良态矩阵
A_good = diag([1, 2, 3, 4]);
eigenvals_good = eig(A_good);

% 添加小扰动
epsilon = 1e-6;
E = epsilon * randn(size(A_good));
A_perturbed = A_good + E;
eigenvals_perturbed = eig(A_perturbed);

fprintf('良态矩阵的特征值扰动:\n');
fprintf('扰动大小: %e\n', norm(E, 'fro'));
eigenval_changes = sort(eigenvals_perturbed) - sort(eigenvals_good);
fprintf('特征值变化: [%.2e, %.2e, %.2e, %.2e]\n', eigenval_changes);

% 病态矩阵（接近奇异）
A_ill = [1, 1; 1, 1+epsilon];
eigenvals_ill = eig(A_ill);

E_ill = epsilon * randn(size(A_ill));
A_ill_perturbed = A_ill + E_ill;
eigenvals_ill_perturbed = eig(A_ill_perturbed);

fprintf('病态矩阵的特征值扰动:\n');
fprintf('条件数: %.2e\n', cond(A_ill));
eigenval_changes_ill = sort(eigenvals_ill_perturbed) - sort(eigenvals_ill);
fprintf('特征值变化: [%.2e, %.2e]\n', eigenval_changes_ill);

%% 9. 应用示例：振动系统
fprintf('\n9. 应用示例：振动系统\n');
fprintf('质量-弹簧系统的固有频率\n\n');

% 质量矩阵和刚度矩阵
M = diag([1, 2, 1]);  % 质量矩阵
K = [2, -1, 0; -1, 2, -1; 0, -1, 1];  % 刚度矩阵

fprintf('质量矩阵 M =\n');
disp(M);
fprintf('刚度矩阵 K =\n');
disp(K);

% 广义特征值问题: K*v = λ*M*v
% 转换为标准特征值问题: M^(-1)*K*v = λ*v
A_vibration = M \ K;
[V_vib, D_vib] = eig(A_vibration);
frequencies = sqrt(diag(D_vib));  % 固有频率

fprintf('固有频率 (rad/s): [%.4f, %.4f, %.4f]\n', frequencies);

% 绘制振型
figure('Name', '振动系统的振型');
for i = 1:3
    subplot(1, 3, i);
    mode_shape = real(V_vib(:, i));
    bar(1:3, mode_shape);
    title(sprintf('第%d阶振型 (f=%.3f)', i, frequencies(i)));
    xlabel('质量点');
    ylabel('振幅');
    grid on;
end

%% 10. 特征值的数值计算挑战
fprintf('\n10. 特征值的数值计算挑战\n');
fprintf('Wilkinson多项式：特征值计算的经典难题\n\n');

% 简化的Wilkinson矩阵
n = 5;
W = diag(1:n) + diag(ones(n-1, 1), 1);
fprintf('Wilkinson型矩阵 (%d×%d):\n', n, n);
disp(W);

eigenvals_W = eig(W);
fprintf('特征值: [%.6f, %.6f, %.6f, %.6f, %.6f]\n', eigenvals_W);

% 条件数分析
cond_W = cond(W);
fprintf('条件数: %.2e\n', cond_W);

% 特征值的条件数（简化估计）
[V_W, D_W] = eig(W);
eigenval_condition_numbers = zeros(n, 1);
for i = 1:n
    v = V_W(:, i);
    w = (V_W' \ eye(n, 1))';  % 左特征向量的近似
    eigenval_condition_numbers(i) = 1 / abs(w * v);
end

fprintf('特征值条件数估计: [%.2e, %.2e, %.2e, %.2e, %.2e]\n', ...
        eigenval_condition_numbers);

end